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Abstract 

In dealing with thermal transport in composite systems, high contrast materials pose a special 
problem for numerical simulation: the time scale or step size in the high conductivity material must 
be much smaller than in the low conductivity material. In the limit that the higher conductivity 
inclusion can be treated as having an infinite conductivity, we show how a standard random walk 
algorithm can be alterred to improve speed while still preserving the second law of thermodynamics. 
We demonstrate the principle in a ID system, and then apply it to 3D composites with spherical 
inclusions. 
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I. INTRODUCTION 



A variety of systems display dramatically different thermal diffusivities. For example, 
the thermal conductivity is estimated at 3000 W/mK for an isolated multiwall carbon nan- 
otube (CNT) and between 1750 and 6600 W/mK for a single wall carbon nanotube at room 

nn 

temperature. [1H4| Typical polymer matrices, in contrast, have thermal conductivities that 
are three orders of magnitude smaller. Composites using carbon nanotubes have been sug- 
gested as cheap materials with average thermal conductivity. [5j However making improved 
thermal conducting polymer composites has been hindered due to the Kapitza thermal re- 
sistance and various processing issues. 

It may be possible to alter this resistance through functionalizing the ends or surface 
of the CNTs, although this may decrease the thermal conductivity of the material. The 
problem of optimizing the thermal conductivity of CNT composites presents an intriguing 
combination of high conductivity contrast, strong disorder, and incorporated materials with 
an extremely high aspect ratio. Thus an efficient and reliable method to calculate effective 
thermal conductivity and varying interface resistance in this two phase medium is desirable. 

This problem has been studied in the context of electrical conductivity. However, 
the inclusions in thermal transport can have be quite asymmetric and entangled. In these 
cases it the isotropic averaging models may not apply. 

One approach to model the thermal transport in composites is to use a random walk 
algorithm in which the transport is assumed to be diffusive. For instance, Tomadakis and 
Sotirchos js, 9| has used this approach to find the effective transpo rt p roperties of random 
arrays of cylinders in a conductive matrix. Recently, Doung et.al 10Hl2| developed a random 
walk algorithm to model thermal transport in carbon nanotube-polymer composites and the 
simulation results showed a reasonable agreement with the experimental data for Epoxy- 



SWNT composites 



111 ]. In this approach thermal transport is described by random jumps of 



thermal markers carrying a certain amount of energy (AE). The step size(Aa; = £2 —^1) of 
this thermal markers follows the gaussian distribution. (See EqJT]) The standard deviation (a) 



of the gaussian step distribution in each one of the space dimensions is a M /j = yj2DM/iAt, 
where At is the time increment, M/I refers to matrix or inclusions and Dm/i is the thermal 
diffusivity. However problem arises when their is a high contrast in thermal diffusivity of 
matrix and inclusions. The step size in highly conducting inclusions Axj become very large 
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compare to the that of poorly conducting matrix Axm- Eventually this leads the markers 
to jump out the inclusions as soon as they enter. This can be avoided having very small 
steps inside the matrix so that the steps inside the inclusions are within the dimensions of 
the inclusions. But this is computationally expensive. 

When the thermal diffusivity of the inclusion is very high relative to the matrix it is 
reasonable to assume that the thermal diffusivity of the inclusions is infinite. This obviates 
the need to model random walks inside the inclusions. In this approach markers entering 
an infinite conductivity inclusion (ICI) are distributed uniformly inside the inclusion on 
the next time step. Some fraction will leave on the next time step and they always leave 
from the surface of the inclusion. (Otherwise the simulation wastes time on walkers that 
hop within the ICI.) However, we must be careful in choosing how the walkers leave the 
the ICI since incorrect approaches can lead to the unphysical result of a system at uniform 
temperature spontaneously developing a temperature gradient at the interface between the 



inclusion and the medium. While the effect is apparently small, it must be remembered 
that diffusion occurs at these same interfaces. In this paper we provide a rigorous approach 
for implementing a random walk algorithm with emphasis on the treatment at the interface 
between the inclusions and the matrix material for high conductivity contrast composites, 
and we quantify the errors made when gaussian and modified step distributions are employed. 

This paper is divided into four parts. In the first, we briefly describe the algorithm for 
"infinite conductivity" inclusions. Next, we show the rigorous way to handle inclusions in one 
dimensional systems. We verify our results numerically in ordered and disordered systems, 
and compare them to results obtained by assuming that the walkers leave the surface with a 
gaussian step distribution. In the next section we develop this approach to spheres in three 
dimensions and again verify it numerically, showing quantitatively the errors that develop 
if a gaussian step distribution from the surface is used. Interestingly, the errors in thermal 
conductivity are larger in 3D than in ID and larger for random arrays than for regular ones. 
In the final section we conclude with a summary and a discussion of future work. 

II. THE MODEL 

The goal is to calculate the thermal conductivity of a composite composed of a matrix 
containing a distribution of "infinite conductors" (ICs). This conductivity is calculated 
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by fixing the heat flux through the computational volume and measuring the resulting 
average temperature gradient. The diffusion of heat is modelled by the motion of random 
walkers within the domain. The computational cell is divided in bins, and the temperature 
distribution is calculated from the number of walkers in each bin. To maintain a constant 
heat flux in the x direction through the computational cell, random walkers carrying +AE 
energy are periodically added at the surface x = x m i n , and then allowed to move with 
random jumps that follow a gaussian distribution into the computational cell. In order 
to fix an "outward" energy flux on the opposite surface, random walkers carrying —AE 
energy are added at the surface x = x max at the same rate as the positive markers. The 
+AE and —AE thermal markers are often called "hot" and "cold" walkers. The exact size 
of AE is arbitrary: the heat flux might be modelled by many small walkers or one large 
one. However using too many walkers is computionally ineffcient, while too few produces 
noisy results that requires more runs to get better averages. In the y and z direction the 
computaional domain is assumed to be periodic. The solution at steady states yields a linear 
temperature profiile and the thermal conductivity can be extracted from Fourier's law. To 
incorporate the effect of the Kapitza thermal resistance, walkers in the matrix that would 



normally attempt to jump into the IC can only do so with a probability f m r^.|13l-ll6| Thus 
they stay in the matrix phase with a probability 1 — f m ,ic- The value of f m ,ic is determined 
by the Kapitza resistance. This can be estimated using acoustic mismatch model when the 
physical properties of the materials are known. 17] 

Similarly, random walkers located within the IC have a probability to hop out on each 
time step. Exactly what fraction of the walkers should leave in each time step, and the 
exact nature of the probability distribution for the steps they should take from the surafce 
are determined in the next two sections. However, those that do leave, exit at random 
positions on the IC. This is done to model the "infinite" conductivity of IC so that the 
walker distribution within the IC is uniform. 

Collisions between walkers are ignored. The random walk reflects the scattering of 
phonons in the disordered matrix material. Walker-walker scattering would reflect nonlinear 
thermal conductivities which are typically small. Similarly, we assume that the properties 
of the materials (e.g. density, specific heat, thermal relaxation length) do not change with 
temperature over the range modelled. 

Finally, we assume that the product of the mass density of IC and specific heat capacity 
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equals that of the matrix, so that in thermal equilibrium the walker density would be uniform 
inside and out of the IC's. This is done for simplicity, so that the local temperature is simply 
proportional to the difference of the average density of hot and cold walkers. Without this 
assumption we would have to alter the probability of walkers entering and leaving the IC's so 
that in thermal equilibrium, the ratio of average walker density inside the IC to that of the 
matrix equals the ratio of their volumetric heat capacities. Only then would the equilibrium 
walker distribution represent a uniform temperature. 

III. RANDOM WALKS WITH INFINITE CONDUCTIVITY INCLUSIONS IN ID. 



Below we describe how to efficiently handle the random walks in a fashion that satisfies 
the second law of thermodynamics. The difficulty lies in properly handling the random 
walkers that jump out from the high conductivity material. To make the explanation clear, 
we first look at the one dimensional case. We subsequently address the three dimensional 
case for spherical inclusions in section IIVI 

A. Analytic results for one dimensional walks 

We consider a set of random walkers moving in a one dimensional ring, half made from 
an "infinitely conducting material" as show in figO We can view this as a one dimensional 
line with boundaries at x — ±1. We know that in equilibrium the density of random walkers 
throughout the whole ring should be uniform. 

Consider a surface located at shown in figfTJb. The flux of random walkers from 

the left through the surface must equal that from the right. This is not a problem for a 
surface located near the center of the interval. However, if 1 > a > 1 — s, the flux of random 
walkers from right in the matrix medium cannot balance those from the left; there are too 
few of them. The solution lies in that the difference must be made up from random walkers 
leaving the "infinite" conductivity material. If they were distributed uniformly throughout 
the infinite conducting material, their flux would maintain the equilibrium. 

We do not wish to model the inside of the IC inclusions because random walkers within 
them move on a much faster time scale than those outside. We assume that a random walker 
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(a) (b) 

FIG. 1: An illustration of the one dimensional model. In (a) the darker material on the left is of 
"infinite" conductivity. The one dimesional system therefore has two boundaries as shown in (b), 
where we assume a « 1. The difficulty is that in equilibrium the net flux through any surface 
(e.g. the dashed line) must be zero. Walkers hopping through the surface from the shaded region 
in Fig.(b) must be balanced by walkers leaving the right hand infinite conductor. This implies 
that if random walkers always leave the surface of the infinite conductor (IC), they must have a 
different jump distribution than random walkers inside the "normal" region. 

instantly leaves from any point on the surface (in this case from x = ±1). However, since 
they leave always exactly from the surface, their step distribution must be different from 
that of random walkers within the matrix medium. 

In each step of the simulation we move the walkers inside the interval — 1 < x < 1, as well 
as those outside. We wish to do this in a fashion that is in agreement with the second law 
of thermodynamics. Let the probability that a walker in the matrix medium jumps from X\ 
to X2 be given by: 

P(x u x 2 ) = -=^e—&— (1) 
V 2nu 

In each time interval we can see that only a fraction of the walkers inside the IC will leave, 
or else their density would not equal that in the normal medium. In this simple model, we 
require that the number inside the matrix medium (Nj) equals that outside (N ) in the IC. 
We also require that the net flux through a surface located at x = s is zero. 

The flux to the left from particles lying in the matrix region, s < x < 1 is balanced by the 
flux to the right for those lying between 2s — 1 < x < s. However those in the shaded region 
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of figO] are not so compensated. They must be balanced by a net flux of walkers leaving the 
righthand boundary. Denote the flux from the shaded region to the right by N r ; it is : 



/2s— 1 roc 
^ dxi J dx 2 P(xi,x 2 ) (2) 
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where Erfc (x) is the complementary error function, Erfc (x) = 1— Erf (x). We are summing 
over any walker starting in the shaded region ending up anywhere to the right of the barrier. 
We have let the upper limit of the endpoint of the jump to infinity since a « 1; we extend 
the lower limit of the first integral to — oo, and we have shifted variables to z = 2s — 1 — x 2 - 
We neglect any walkers leaping from the IC on the left boundary x — — 1 all the way through 
s. 

The flux N r must be balanced by the flux from walkers leaving the IC on the right. Let 
the probability that a random walker in the IC leaves it be given by A, and the probability 
that it jumps to a point x, leaving from the right hand boundary, be f{x). Then the flux 
to the left through the surface at x = s due to these walkers is 

N e = N X f f(x)dx (4) 

J — oo 

We set Ne = N r , and take the derivative of both sides with respect to s. This gives us an 
integral expression for f(s): 



m p0 



2N X 



(5) 



The requirements that Ni = N and the balancing of the fluxes when s = 1 is enough to 
solve for f(s). The distribution of steps, f{u) = f{l — u) is given by: 

B. Numerical results for thermal conductivity in ID 

The above analytical calculation provides the correct step distribution for walkers leaving 
the edge of the infinite conductors. We can compare it to a simple model where we simply 
have the walkers take a step with a Gaussian probability distribution (mean size 0.20) from 
the surface. Figj2] is the spatial distribution of random walkers in such a one dimensional 
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FIG. 2: (Color online) Results for the one dimensional model of a random walk on a ring with an 
infinite conductivity inclusion. The steps in the random walk have a Gaussian distribution with 
a mean of 0.20. Plotted are the average number of walkers in each of 20 bins, equally spaced 

— 1 < x < 1 over the course of 10 5 Monte Carlo steps. Starting with 500 random walkers, half 
should be in the "normal" region, so that the average number/bin should be 12.5. The dashed line 
is the result of performing the simulation incorrectly, and letting the walkers have a Gaussian step 
distribution as in eq{TJ There are too many walkers in the interval, and their distribution is not 
uniform. The solid line is the result of using eqlH which yields the correct result. 

system. [3] Plotted are the average number of walkers in each of 20 bins, equally spaced 

— 1 < x < 1 over the course of 10 5 Monte Carlo steps. Starting with 500 random walkers, 
half should be in the "matrix" region, so that the average number/bin should be 12.5. The 
dashed line is the result of performing the simulation incorrectly, and letting the walkers 
have a Gaussian step distribution as in eq{TJ There are too many walkers in the interval, 
and their distribution is not uniform. The solid line is the result of using eqj6l which yields 
the correct result, and is uniform. 

In order to determine the significance of this error, we place several ICs in the compu- 
tational volume and run at constant heat flux until the temperature distribution converges. 
We then extract the gradient in walker density and calculate the thermal conductivity. Sam- 
ple results are plotted in fig.(j3]), where we show the results for Gaussian steps (lower curve) 
and steps governed by eq.@ (upper curve). The latter gives physically reasonable results 
(with noise), in which the temperature is constant the ICs and uniformly decreasing in the 
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FIG. 3: (Color online) Plots of the temperature distribution in a periodic array of infinite conduc- 
tivity inclusions in a ID matrix. The temperature is determined by the average of the difference 
between the number of positive and negative walkers in a given region. This plot was generated 
over a run of 62,500 time steps where 10,000 positive (negative) walkers were added every 10 time 
steps at the left (right) border. The ICs are 0.50 units long and the matrix spacer between them is 
1.00 units. The probability to cross into an IC from the matrix is 0.16. Lower line: The tempera- 
ture distribution generated by having random walkers depart the ICs with Gaussian steps. Upper 
line: The temperature distribution (shifted up one unit for clarity) given by an algorithm using 
the step probability distribution of eqEJ Note that temperature gradients spontaneously appear 
at interfaces when the incorrect jump distribution is used. 

matrix. 

The thermal conductivity is extracted from the ratio of the slope of the temperature (the 
walker density) to the applied flux. In fig. (jH ) we plot the average value of the percent error 
in the thermal conductivity as a function of the transmission probability, f m jc, for regular 
and random ID arrays. In this simulation the ICs were 0.50 units long and the material 
between them was 1.00 units wide. The results are averaged over five runs each lasting 
for 40,000 time steps. The percent error is denned as the difference between the results of 
simulations using the Gaussian steps and the results using eqj6l The error bars represent 
the variation in thermal conductivities over the runs. Thus we see that the error can range 
as large as five percent, and that can vary substantially. 
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FIG. 4: (Color online) Plot of the relative percent error for the the thermal conductivity calculated 
using a Gaussian step distribution (a = 0.1) as compared to that of eq.© as a function of f m -ic, 
the probability for a walker to enter into an inclusion. The results are for a one dimensional system 
with 20 inclusions; the diamonds are for a regular array of ICs and the circles are random ICs. 
The error bars are based on a sample of five different configurations and are included to give an 
indication of how large the the errors can be. 

IV. THREE DIMENSIONAL MODELS 

We have shown above that errors in one dimensional simulations are avoidable, but only 
a few percent. Below we generalize the above problem to three dimensions for spherical 
inclusions, and show that the effect can be significant. 

A. Analytic Derivation of Random Walks with Spherical Inclusions 

In our model random walkers that land inside the sphere are immediately moved to a 
random point on the surface of the sphere. On the next time step they can move in the radial 
direction away from the sphere. We assume that if we choose the fraction that leave and 
their step distribution correctly, then when we are in equilibrium we will obtain a uniform, 
stationary density outside and inside the sphere. 

The number of random walkers entering the sphere from a region dr near r landing inside 
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the sphere is given by: 



n{r) = dr 



Pofk 



-(r-r') z 



dr' e 



(7) 



(2tt) 3 /V 3 Jr><R 

where we have generalized the probability distribution of eq.((T]) to three dimensions. The 
factor fk represents the Kapitsa resistance; it is the probability that a random walker will 
enter the spherical inclusion. When /& = 0, no walkers enter the inclusion and the Kapitsa 
resistance is infinite; when f% = 1 then walkers can freely step into the inclusion and the 
Kapitsa resistance is zero. The total number entering a sphere of radius R is 

Pofk 



N in (R) 



dr 



-(r-r'r 

dr' e z? 2 



(8) 



r>R (27r) 3 / 2 <7 3 Jr'<R 

For any value of r we can rotate our primed coordinate system so that z! || f so that the 
angle between r and r ' is simply 9', the spherical polar angle in the primed system. The 
angular integrals can then all be done in closed form giving 
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r dr 
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r 'dr' 
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e 2^ 



(9) 



The resulting integral can also be found exactly yielding 
2 



N in (R) = vPofk 



2ttct 3iT - a' 



2ttR 3 Erfc (^-^) + ^phxae 



2R' 



(a 2 - R 2 



(10) 



If the density of walkers is uniform then the number inside the sphere is V s po where V s 
is the volume of the sphere. (This is only true when the product of the density and specific 
heat capacity of the matrix and IC's are equal. When this constraint does not hold the 
number of walkers inside the IC is V s p c^pw ' ^ nere Cic(Cm) and Pw(pm) are specific 
heat capacity and mass density of the IC (Matrix).) In each time step we allow a fraction A 
of them to leave. In equilibrium the flux into the sphere (Eq fl0|) equals the flux out (V s po\), 
allowing us to calculate A: 



A 



111 
3K 



27r<r(3ir - a 2 ) + 27riT erfc 



a 



2gf 

) + V2nae ^ (cr 



R 2 



(11) 



When R » a, we expect the geometry of the inclusion to be irrelevant. In this limit, if 
the random walkers had a flat distribution of steps bounded by a, then the flux into the 
sphere would come from a thin spherical shell of thickness a and radius R. The volume of 
this shell is crA s , where A s is the surface area of the sphere. The flow in from this shell is 
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FIG. 5: An illustration of the three dimensional model for spheres. The grey material is a sphere 
of "infinite" conductivity with radius R. We wish to calculate the distribution of steps sizes taken 
by random walkers leaving the surface of the sphere. We do this by requiring that in equilibrium 
the net flux through any surface (e.g. the dashed line indicating a sphere of radius s) must be 
zero. The number of walkers hopping in through the dotted surface (iVj n (s)) must be balanced by 
the total number of walkers hopping out, both those from the matrix material (N ou t(R,s)) and 
those from the surface of the sphere (N sp ^ ere (s)) . This condition allows us to calculate the step 
distribution. 

balanced by the flow out of the volume, V s pX. It is useful to write this in terms of a new 
constant, cq, defined via 

« = % < 12 > 

which is dimensionless and becomes shape independent as a — > 0. In this case 



This quantity is bounded by fk/v2n, the result one would get for an infinite slab. The 

factor l/y/2% arises from the fact that walkers have a gaussian distribution of step sizes, 
and not a flat one. 

Next we have to calculate the distribution of steps for random walkers leaving the surface 
of the sphere. As in the one dimensional case of subsection IIII Al above, we can calculate the 
desired result by balancing fluxes in equilibrium. We draw an imaginary surface of radius 
s about the spherical inclusion. In equilibrium, the net flux through this surface must be 
zero, as illustrated in fig.©. 




(13) 
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N in (s) = N out (R, s) + N sphere (R, s) (14) 

The flux in through the sphere of radius s, N in (s), is the result eq.( [10]) evaluated for a 
radius of s. The flux outward from the matrix material is given by the integral: 

Nout(R,s)=f dr J°J" - [ dr'^ 1 (15) 



lR<r<s (2ir) 3 / 2 a 3 Jr'> 
We can again evaluate this integral analytically to obtain: 
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N ont (R, s) = -f k p 



27rs 3 erfc (— ) + V2^a(a 2 - s 2 )e^ - V2^a(a 2 - 3s 2 ) 
a 



+V27ia((r z — s — R — Rs)e^^ - v27rcr(<7 — s — R + Rs)e ^ 
+7r(s 3 - R 3 ) erfc f^5) - vr( S 3 + R 3 ) erfc r S R 



(16) 



\/2a V2cj 

Finally, we can write an expression for the flux of random walkers (originating on the 
inclusion surface) that hop out through the sphere of radius s: 

f(r)dr (17) 

where f(r) dr gives the fraction of walkers that jump radially outward to a distance between 
r and r + dr from the center of the inclusion. 

Eqns.(|T0|). (|T6|) and (|T7|) give us enough information to calculate the step distribution 
function f(r). However in computer applications we do not actually use f(r). Rather 
algorithms typically generates a random number, p, in a flat distribution < p < 1, and use 
that to select a random step S(p) from the center of the sphere, 5 > R. We can do this by 
first calculating the integral of f(r): 

P(8) = f f(r') dr' (18) 

J Ft 

Note that P{R) = and lim^oo P(5) = 1. We then must invert this functional relationship 
to get S(P), which gives us the step generating function we desire. We note that from eq. (|T7|) 
and (|T8|) we have: 

N sphcre (R, s) = V s Po A [1 - P(s)} (19) 



Equating this via eq. (1141) and dropping exponentially small terms we have 
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P(s) = 1- 



n(R 3 - s 3 ) erfc [ S -^S\ - V^a(a 2 - s 2 - R 2 - Rs)e 



V2a 

+ v 7 ^^ 2 -s 2 -R 2 + Rs)e^ }1 + n(s 3 + R 3 ) erfc (^=g) 



v / 2^a(3i? 2 - a 2 ) + v / 2^a(a 2 - R 2 )e^~ + 27T.R 3 erfc 



fy/2R\ 



(20) 



V a / 

This result has the desired behavior at the limits, P{R) = and lim^oo P(s) = 1. This 
function is not analytically invertible; in implementation it is evaluated on a mesh and the 
inverse is calculated via interpolation. 

B. Numerical Results in 3D 

We implemented a random walk algorithm in three dimensions similar to that of section 
IIHI above. In three dimensions we applied periodic boundary conditions in the y and z 
directions. A temperature profile in the x direction was obtained simply by binning all 
walkers in a given range of x for all y and z\ such slices would cross inclusions as well as 
matrix material. Walkers that were labelled as inside a given inclusion were assigned a 
random position inside the inclusion for the purpose of doing this averaging. The simulation 
volume was 10 x 10 x 10, and the random walk in the matrix was described by a Gaussian 
distribution with a rms value of 0.10 in these units. The transition probability f m ,ic was 
fixed at 1.0. 

In fig. ([6]) the percent error (defined as the ratio of the difference of thermal conductivities 
measured using the Gaussian step distribution and that of eq.f l20|) . divided by the former) 
is plotted as a function of the volume fraction of infinite conductivity inclusions, for a fixed 
surface area for the inclusions. (If there were only a single spherical inclusion, it would 



have had a volume fraction of 5%.) 19] As the number of inclusions at fixed surface area 
increases, their total volume decreases as N~ 3 / 2 . (For example, the largest volume fraction, 
0.20, corresponds to 100 spheres of radius 0.9772). The results at several values of iV were 
calculated for five random configurations and the average and standard deviation are plotted. 
The effect of using a simplified step distribution is larger in three dimensions, and can affect 
the results by up to 18%. 
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FIG. 6: Plot of the percent error in the thermal conductivity as a function of the volume fraction for 
fixed surface area. The percent error is defined as the ratio of the difference of thermal conductivities 
measured using the Gaussian step distribution and that of eq. (|20p divided by the former. Note 
that the error varies only slightly with volume fraction. 

The percent error was also calculated as a function of the surface area for fixed volume 
fraction and plotted in fig. ([7]). The volume fraction was fixed at 5%, and the surface area 
increases with N as N 2 ^ 3 . Five simulations were run for N = 100, 200, ... 1000 and the 
average and standard deviation were plotted as a function of the surface area relative the 
minimum surface area, A , the area of a sphere that is 5% of the volume. Again, the effect 
of using the wrong simulation algorithm is shown to be substantial. 

V. CONCLUSIONS AND FUTURE WORK 

Transport in composites with a large disparity in conductivities is important to a large 
number of systems. In this paper we have demonstrated an efficient and physically sound 
algorithm for calculating effective conductivities of composites with large contrasts in con- 
ductivity. We have shown that the errors introduced are small but measurable in one 
dimension, and moderately significant in 3D. 

The spherical inclusion case is the simplest 3D problem, but not the most relevant to many 
systems. Carbon nanotubes might be approximated as cylinders, to lowest order. However, 
in that case the ID integrals of section (1IV Aj) become more complicated and handling the 
endcaps of the cylinders becomes problematic. A simple approach might be to simply ignore 
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FIG. 7: Plot of the percent error in the thermal conductivity as a function of the surface area of 
the spherical inclusions at a fixed volume fraction of 5%. The surface area is measured in terms 
of Aq, the surface area of a sphere with 5% of the total volume. The percent error is defined as 
the ratio of the difference of thermal conductivities measured using the Gaussian step distribution 
and that of eq. (j20j) . divided by the former. 

transport through the endcaps, or treat a nanotube as an extremely prolate spheroid so that 
diffusion from the inclusion can again be treated as a one dimensional walk normal to the 
surface. These approximations are the subject of current research. 
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